% File: sc_fde_base.m
% -------------------
% This script is the simulation main body.

% -------------------------------------------------------------------------
% test point - mse, var buffer
% -------------------------------------------------------------------------
mse_buffer      = zeros(length(SNR),5,nsubf); % buffering mse calculate, 5 iterations
mse_test_buffer = zeros(size(mse_buffer));
var_buffer      = zeros(size(mse_buffer));
var_test_buffer = zeros(size(mse_buffer));
% -------------------------------------------------------------------------

numerr    = zeros(length(SNR),nsubf); % Compute error rate for each subframe
bitnumerr = zeros(size(numerr));
ser_mfb   = zeros(1,length(SNR));
ber_mfb   = zeros(size(ser_mfb));
h_waitbar = waitbar(0,'Please wait...');
for n = 1:length(SNR)
%     chan_mimo = mimochan_init(nt,nr,Ts,doppler_max,path_delays_norm,pdp);
    mimochan_init          % mimo channel initialization, refers ot mimochan_init.m for details
    m = 1;                 % reset loop counter
    while (m <= numTrials) % quasi-static channel modeling
        switch chan_profile
            case 1  % fixed channel for text only
                h = h_test_rayleigh;
                h_sq_norm = sum(abs(h).^2); % channel pulse response energy
%                 H = fft(h_test_rayleigh, nfft, 2).';
            otherwise
                ... no op.
        end
        %%
        ind_perm_real = zeros(nfft,nsubf);                % chip-wise interleaving index for real part
        ind_perm_imag = zeros(nfft,nsubf);                % chip-wise interleaving index for imag part
        for p = 1: nsubf
            if(on_off_interleaving_chip)
                ind_perm_real(:,p) = randperm(nfft).'; % random interleaving
                ind_perm_imag(:,p) = randperm(nfft).'; % random interleaving
            else
                ind_perm_real(:,p) = (1:nfft).';       % no interleaving
                ind_perm_imag(:,p) = (1:nfft).';       % no interleaving
            end
        end
        if(on_off_freq_interleaving_chip)              % chip-wise frequency interleaving index
            ind_chip_perm_freq = randperm(nfft).';     % random interleaving
        else
            ind_chip_perm_freq = (1:nfft).';           % no interleaving
        end
        %%
        if (~code_profile)
            x     = randi([0,M_mod-1], nsymb, nsubf); % create a random digital message, use M-QAM modulation.
        else
            xbit  = randi([0 1], nbit, nsubf);
            xcode = zeros(nbit/crate, nsubf);
            x     = zeros(nsymb, nsubf);
            ind_bit_perm = zeros(size(xcode));
            for p = 1: nsubf
                xcode(:,p) = convenc(xbit(:,p),trel);               % encoding
                ind_bit_perm(:,p) = randperm(size(xcode,1));        % bit-wise interleaving index
                xcode(:,p) = intrlv(xcode(:,p), ind_bit_perm(:,p)); % interleaving
                x(:,p)     = bi2de(reshape(xcode(:,p), bit_sym, length(xcode(:,p))/bit_sym).','left-msb'); % bit-to-symbol mapping
            end
        end
        if (on_off_scrambling) % scrambling, no need for de-scrambling at rx because the scrambled symbols are the new signal in concern
            calc_sc = M_mod;   % calculation base
            poly_sc = [6 5 0]; % scramble polynomial
            init_sc = randi([0,M_mod-1], 1, max(poly_sc)); % initial states of scrambler's registers
            for p = 1: nsubf
                x(:,p) = scrambler(x(:,p), calc_sc, poly_sc, init_sc);
            end
        end
        x_pload = modulate(modem.qammod(M_mod),x);
        x_pload_spread = repmat(x_pload,[sf,1]); % symbol spreading into chips and interleaving
        xtx_pload      = zeros(size(x_pload_spread));
        for p = 1: nsubf
            xtx_pload(:,p) = real(x_pload_spread(ind_perm_real(:,p),p))...
                             + 1i*imag(x_pload_spread(ind_perm_imag(:,p),p));
        end
        XTX   = fft(xtx_pload,nfft,1);     % DFT of transmitted symbols, used at rx end for filter design.
        XTX_p = XTX(ind_chip_perm_freq,:); % chip-wise frequency interleaving
        xtx_p = ifft(XTX_p,nfft,1);
        %%
        if (length_cp>0)
            xtx_subf = [xtx_p(end-length_cp+1: end, :); xtx_p]; % insert CP
        else
            xtx_subf = xtx_p; % no CP at all
        end
%        xtx = [xtx_pload; zeros(length_cp,1)]; % add UW at the end of a block.
%                                          % the UW at the beginning of a
%                                          % block is omitted, since the UW
%                                          % we choosed here is a null
%                                          % sequence, it will not interfere
%                                          % with the pay load.
        %%
        xtx = reshape(xtx_subf,[],1);  % superimpose sub-frames into a superframe
        switch chan_profile
            case 1
                fadesig = filter(h,1,xtx);
                H = fft(h, nfft);
                H(1,ind_chip_perm_freq) = H;
            otherwise
                %[fadesig h h_sq_norm] = mimochan_filter(chan_mimo, xtx); 
                mimochan_filter                % mimo channel filtering, refers to mimochan_filter.m for details
                H   = fft(h, nfft, 3);         % channel estimation, ideal
                H(:,:,ind_chip_perm_freq) = H; % enable for chip-wise frequency interleaving
        end
        %%
        xtxPower      = sum(abs(xtx).^2,1)/size(xtx,1);   % check transmitted signal power
        xtxPower_mean = mean(xtxPower);
        xtxPower_mean_dB = 10*log10(xtxPower_mean);       % convert to dB
        noisePower_dB = xtxPower_mean_dB-SNR(n);          % compute noise power
        noisePower    = 10^(noisePower_dB/10);            % convert backto linear scale
        noisePow_norm = noisePower/xtxPower_mean;         % noise power normalized to tx power
        rnoise = awgn(fadesig, SNR(n), xtxPower_mean_dB); % additive noise with the assumption
                                                          % that the multipath channel gains
                                                          % are normalized
        % decompose received superframe into sub-frames
        rnoise_subf = reshape(rnoise,[],nsubf);
        %%
        if (length_cp>0)
            rnoise_subf = rnoise_subf(length_cp+1:end, :); % remove CP
                                                           % for UW case, no need to remove anything
        end
        fdein = fft(rnoise_subf,nfft,1);     % convert into frequency domain for FDE
        fdein(ind_chip_perm_freq,:) = fdein; % chip-wise frequency de-interleaving to be inserted here
        %%
        switch rx_profile
            case 1
                sc_fde_ib               % baseline iterative equalizer
            case 2
                sc_fde_rx_cprecon
            case 3
                sc_fde_ib_maxvar
            case 4
                sc_fde_ib_geomvar
            case 5
                sc_fde_ib_harmvar
            case 6
                sc_fde_ib_qam           % iterative equalizer for qam constellation, simple estimate
            case 7
                sc_fde_ib_qam_evo       % iterative equalizer for qam constellation, evolutional estimate
            case 8
                sc_fde_ib_qam_genie     % iterative equalizer for qam constellation, genie aided estimate
            otherwise
                z = x;                       % dummy rx, mfb curve only
        end
        %%
        % Performance evaluation
        if (isempty(find(isnan(z(:)),1))) % check if equalizer outputs NaN values
            [num]       = symerr(z,x,'column-wise'); % symbol error rate in each frame
            numerr(n,:) = numerr(n,:)+num;
            if (code_profile)
                zcode = zeros(size(xcode));
                zbit  = zeros(size(xbit));
                for p = 1: nsubf                     % symbol-to-bit mapping
                    zcode(:,p) = reshape(de2bi(z(:,p),'left-msb').',[],1);
                    zcode(:,p) = deintrlv(zcode(:,p), ind_bit_perm(:,p));          % de-interleaving
                    zbit(:,p)  = vitdec(zcode(:,p), trel, tbleng, 'cont', 'hard'); % Decode.
                end
                [bitnum]       = biterr(zbit(decdel+1:end,:), xbit(1:end-decdel,:), 'column-wise');
                bitnumerr(n,:) = bitnumerr(n,:)+bitnum;
            end
            
            Cr_mfb  = 1; % Compute MFB
            sf_mfb  = 1;
            SNR_MFB = Eb_N0(n) + 10*log10(Cr_mfb*log2(M_mod)) - 10*log10(sf_mfb);
            snr_mfb = 10^(SNR_MFB/10)*h_sq_norm;
            kap = 1; % kap = (.5*d_min)^2/((i^2+q^2)/N); where d_min is the 
                     % minimum distance between a constellation pair, and N
                     % equals to 2 when using quadrature (IQ) modulation.
            mfb = snr_mfb*kap;
            ser_mfb(n) = ser_mfb(n) + 2*.5*erfc(sqrt(mfb)/sqrt(2)); % Q(x) = .5*erfc(x/sqrt(2))
            ber_mfb(n) = ber_mfb(n) + .5*erfc(sqrt(mfb)/sqrt(2)); % Q(x) = .5*erfc(x/sqrt(2))
            
            waitbar(((n-1)*numTrials+m)/(numTrials*length(SNR)))
            
            m = m+1;
        else
            ... % keep m unchanged if equalizer outputs NaN values
        end
    end
end
close(h_waitbar)

% End of script